module dg_tools
implicit none
    real, parameter, dimension(4),private :: gauss_points = [ &
        0.3399810435848562648026658, &
        0.8611363115940525752239465, &
        -0.3399810435848562648026658, &
        -0.8611363115940525752239465 ]
    
    real, parameter, dimension(4),private :: gauss_coef = [ &
        0.6521451548625461426269361, &
        0.3478548451374538573730639, &  
        0.6521451548625461426269361, &
        0.3478548451374538573730639 ]   
    
    !pi值
    real, parameter :: pi = 3.141592653589793238462643
contains
    ! 输出高斯积分点和权值
    subroutine show_consts()
        print '(A)', "Gauss_coefficient:"
        print '(F25.20)', gauss_coef  ! 高斯积分权值
        print '(A)', "Gauss_points :"
        print '(F25.20)', gauss_points ! 积分点
    end subroutine show_consts
    
    ! 获取高斯积分点
    subroutine get_points(points)
        implicit none
        real, intent(out), dimension(4) :: points
        points = gauss_points
    end subroutine get_points
    subroutine get_pi(ppi)
        real, intent(out) :: ppi
        ppi = pi
    end subroutine get_pi
    ! 计算基函数在高斯积分点处的值
    subroutine compute_bvalue(k, b_value)
        implicit none
        integer, intent(in) :: k
        real, intent(out), dimension(k+1,4) :: b_value
        integer :: n
        select case (k)
            case (1)
                do n = 1,4
                    b_value(1,n) = 1
                    b_value(2,n) = gauss_points(n)
                end do
            case (2)
                do n = 1,4
                    b_value(1,n) = 1
                    b_value(2,n) = gauss_points(n)
                    b_value(3,n) = 0.5*(3*(gauss_points(n)**2) - 1)
                end do
        end select
    end subroutine compute_bvalue
    subroutine compute_dbvalue(k, db_value)
        implicit none
        integer, intent(in) :: k
        real, intent(out), dimension(k+1,4) :: db_value
        integer :: n
        select case (k)
            case (1)
                do n = 1,4
                    db_value(1,n) = 0
                    db_value(2,n) = 1
                end do
            case (2)
                do n = 1,4
                    db_value(1,n) = 0
                    db_value(2,n) = 1
                    db_value(3,n) = 3*gauss_points(n)
                end do
        end select
    end subroutine compute_dbvalue
    ! 高斯积分
    subroutine compute_int(fun_value,int_value)
        real, intent(in), dimension(4) :: fun_value
        real, intent(out) :: int_value
        int_value = dot_product(gauss_coef, fun_value)
    end subroutine compute_int

    !计算基函数在边界点以及中点的值
    subroutine compute_basis(k,basis)
        integer,intent(in) :: k
        real,intent(out), dimension(k+1,3) :: basis
        basis(1,:) = 1.0
        basis(2,1) = -1.0
        basis(2,2) = 1.0
        basis(2,3) = 0.0
        if(k==2) then
            basis(3,:) = 1.0
        end if
    end subroutine compute_basis
    subroutine compute_dbasis(k,dbasis)
        integer,intent(in) :: k
        real,intent(out), dimension(k+1,3) :: dbasis
        dbasis(1,:) = 0.0
        dbasis(2,:) = 1.0
        if(k==2) then
            dbasis(3,1) = -3.0
            dbasis(3,2) = 3.0
            dbasis(3,3) = 0.0
        end if
    end subroutine compute_dbasis
    !生成质量矩阵
    subroutine mass_matrix(k,h,Mass)
        real, intent(in) :: h
        integer, intent(in) :: k
        real, intent(out) ,dimension(k+1,k+1) :: Mass
        real :: c
        c = 2.0 / h
        Mass(1,1) = c*0.5
        Mass(2,1) = 0.0
        Mass(1,2) = 0.0
        Mass(2,2) = c*1.5
        if (k == 2) then
            Mass(1,3) = 0.0
            Mass(2,3) = 0.0
            Mass(3,3) = c*2.5
            Mass(3,1) = 0.0
            Mass(3,2) = 0.0
        end if
    end subroutine mass_matrix
    ! 生成刚度矩阵
    subroutine stiff_matrix(k,stiff)
        integer, intent(in) :: k
        real, intent(out), dimension(k+1,k+1) :: stiff
        stiff(:,1) = 0.0
        stiff(:,2) = 0.0
        stiff(2,1) = 2.0
        if (k==2) then
            stiff(:,3) = 0.0
            stiff(3,2) = 2.0
        end if
    end subroutine stiff_matrix
    !初值投影
    subroutine get_initial_coef(k,cell_num,h,u_value,a_coef,mass)
        integer,intent(in) :: k, cell_num
        real,intent(in) :: h
        real,intent(in), dimension(cell_num,4) :: u_value
        real,intent(out), dimension(k+1, cell_num) :: a_coef
        real,intent(in), dimension(k+1,k+1) :: mass
        real, dimension(k+1,4) :: b_value
        real, dimension(4) :: f_value
        real, dimension(k+1) :: right_b
        integer :: n,i
        call compute_bvalue(k, b_value)
        do n =1 , cell_num
            do i =1 ,k+1
                call vector_mul(4,b_value(i,:),u_value(n,:), f_value)
                call compute_int(f_value, right_b(i))
                right_b(i) = right_b(i) * (h/2.0)
            end do
            a_coef(:,n) = matmul(mass, right_b)
        end do
    end subroutine get_initial_coef
    ! 向量对应元素相乘
    subroutine vector_mul(n,a,b,vec)
        integer, intent(in) :: n
        real, intent(in), dimension(n) :: a
        real ,intent(in), dimension(n) :: b
        real, intent(out), dimension(n) ::vec
        integer :: i
        do i =1 ,n
            vec(i) = a(i) * b(i)
        end do
    end subroutine vector_mul
    ! 计算L2误差
    subroutine compute_l2error(h,k,cell_num,a_coef,u_value,error)
        integer,intent(in) ::k, cell_num
        real,intent(in) :: h
        real,intent(out) :: error
        real,intent(in), dimension(k+1, cell_num) :: a_coef
        real,intent(in), dimension(cell_num,4) :: u_value
        real,dimension(k+1,4) :: b_value
        real,dimension(cell_num,k+1) :: a_trans
        real,dimension(cell_num,4) :: num_value
        real,dimension(4) :: tmp
        real,dimension(4) :: error_square
        real :: tmp_int
        integer :: i
        call compute_bvalue(k, b_value)
        a_trans = transpose(a_coef)
        num_value = matmul(a_trans, b_value)
        do i=1, cell_num
            tmp = num_value(i,:) - u_value(i,:)
            call vec_square(4, tmp, error_square)
            call compute_int(error_square,tmp_int)
            tmp_int = (h/2.0)*tmp_int
            error = error + tmp_int
        end do
        error = sqrt(error)
    end subroutine compute_l2error
    ! 计算L1误差
    subroutine compute_l1error(h,k,cell_num,a_coef,u_value,error)
        integer,intent(in) ::k, cell_num
        real,intent(in) :: h
        real,intent(out) :: error
        real,intent(in), dimension(k+1, cell_num) :: a_coef
        real,intent(in), dimension(cell_num,4) :: u_value
        real,dimension(k+1,4) :: b_value
        real,dimension(cell_num,k+1) :: a_trans
        real,dimension(cell_num,4) :: num_value
        real,dimension(4) :: tmp
        real :: tmp_int
        integer :: i
        call compute_bvalue(k, b_value)
        a_trans = transpose(a_coef)
        num_value = matmul(a_trans, b_value)
        do i=1, cell_num
            tmp = abs(num_value(i,:) - u_value(i,:))
            call compute_int(tmp,tmp_int)
            tmp_int = (h/2.0)*tmp_int
            error = error + tmp_int
        end do
        error = sqrt(error)
    end subroutine compute_l1error
    ! 计算阶数
    subroutine compute_order(n,error)
        integer, intent(in) :: n
        real, intent(in), dimension(n) :: error
        integer :: i
        real :: order
        order = 0.0
        do i = 1, n-1
            if (error(i+1) /= 0.0) then
                order = log(error(i)/error(i+1)) /log(2.0)
            end if
            print '(A,F12.8)', "Order of convergence: ", order
        end do
    end subroutine compute_order
    subroutine vec_square(n, vec, result)
        integer, intent(in) :: n
        real, intent(in), dimension(n) :: vec
        real, intent(out), dimension(n) :: result
        integer :: i
        do i = 1, n
            result(i) = vec(i)**2
        end do
    end subroutine vec_square
    subroutine basicCL_one_step_Euler(k, cell_num, a_coef, dt,basis,mass, stiff)
        integer,intent(in) :: k, cell_num
        integer :: nn, l, p
        real, intent(in) ::dt
        real, intent(in), dimension(k+1,k+1) :: mass, stiff
        real, intent(in), dimension(k+1,3) :: basis
        real, intent(inout), dimension(k+1,cell_num) :: a_coef
        real, dimension(cell_num,2) :: u_bound
        real, dimension(k+1) :: b

        u_bound = matmul(transpose(a_coef),basis(:,1:2))
        do nn = 1, cell_num
            if (nn ==1 ) then
                p = cell_num
            else
                p = nn - 1
            end if
            do l = 1, k+1
                b(l) = -u_bound(nn,2)*basis(l,2) + u_bound(p,2)*basis(l,1)
            end do
            a_coef(:,nn) = a_coef(:,nn) + dt * matmul(mass, b) + dt* matmul(mass, matmul(stiff, a_coef(:,nn)))
        end do
    end subroutine basicCL_one_step_Euler
    ! Rk3
    subroutine basicCL_RK3(k, cell_num, a_coef, dt, basis, mass, stiff)
        integer,intent(in) :: k, cell_num
        real, intent(in) :: dt
        real, intent(in), dimension(k+1,k+1) :: mass, stiff
        real, intent(in), dimension(k+1,3) :: basis
        real, intent(inout), dimension(k+1,cell_num) :: a_coef
        real, dimension(k+1,cell_num) :: a_coef1, a_coef2, recover
        recover = a_coef
        call basicCL_one_step_Euler(k, cell_num, a_coef, dt, basis, mass, stiff)
        a_coef1 = a_coef
        call basicCL_one_step_Euler(k, cell_num, a_coef1, dt, basis, mass, stiff)
        a_coef2 = 0.75*recover + 0.25*a_coef1
        call basicCL_one_step_Euler(k, cell_num, a_coef2, dt, basis, mass, stiff)
        a_coef = (1.0/3.0) * recover + (2.0/3.0) * a_coef2
    end subroutine basicCL_RK3
end module dg_tools
program DG_cl
    use dg_tools
    implicit none
    real, dimension(4) :: points
    integer :: k, cell_num,i,j,n
    real :: h,a,x_left,x_right,x_mid,tend,dt,t
    real, dimension(5) :: error
    real, dimension(:,:), allocatable :: u_ini_value
    real, dimension(:,:), allocatable :: a_coef
    real, dimension(:,:), allocatable :: u_value
    real, dimension(:,:), allocatable :: mass
    real, dimension(:,:), allocatable :: stiff
    real, dimension(:,:), allocatable :: basis
    call get_points(points)
    a = 0.0
    k = 2
    tend = 1.0
    allocate(stiff(k+1,k+1))
    allocate(basis(k+1,3))
    call stiff_matrix(k, stiff)
    call compute_basis(k, basis)
    do n = 1,5
        t = 0.0   
        h = 1.0 / (8*(2**n))
        dt = 0.2*h
        cell_num = 8*(2**n)
        allocate(u_ini_value(cell_num,4))
        allocate(u_value(cell_num,4))
        allocate(a_coef(k+1,cell_num))
        allocate(mass(k+1,k+1))
        call mass_matrix(k, h, mass)
        ! 逐cell遍历
        do i =1, cell_num
            x_left = a + (i-1)*h
            x_right = a + i*h
            x_mid = (x_left + x_right) / 2.0
            do j = 1,4
                u_ini_value(i,j) = u_ini(x_mid + 0.5*h*points(j))
                u_value(i,j) = u_ini(x_mid + 0.5*h*points(j)-tend)
            end do
        end do
        call get_initial_coef(k, cell_num, h ,u_ini_value, a_coef,mass)
        ! 推进时间层
        do while(t < tend)
            if (t + dt >tend) then
                dt = tend -t
            end if
            t = t + dt
            call basicCL_RK3(k, cell_num, a_coef, dt, basis, mass, stiff)
        end do
        call compute_l2error(h, k, cell_num, a_coef,u_value, error(n))
        deallocate(u_ini_value)
        deallocate(a_coef)
        deallocate(mass)
        deallocate(u_value)
    end do
    call compute_order(5,error)
    deallocate(stiff)
    deallocate(basis)
    contains
    real function u_ini(x)
        implicit none
        real, intent(in) :: x
        u_ini = sin(2.0*pi*x)
    end function u_ini
end program DG_cl